#replication code for Umeda, Michio "The Politics of Aging: Age Difference in Welfare  Issue Salience in Japan 1972-2016"

data <- read.table("replicationdata_PoliticsOfAgingUmeda_APR2020.txt")

#install.packages("lme4")
library(lme4)

#Main analyses 
model1_welfare <- glm(welfare ~ factor(age5),
                       family = binomial, data = data) 

model2_welfare <- glmer(welfare ~ factor(age5) + year20 + n_article
  + factor(gender) + factor(edu3) + factor(citysize) 
  + factor(issue)+ LDPvote + unemployment
  + factor(wording) + factor(chamber) + (1|year),
  family = binomial, data = data) 

model3_welfare <- glmer(welfare ~ factor(age5)*year20 + n_article
  + factor(gender) + factor(edu3) + factor(citysize) 
  + factor(issue)+ LDPvote + unemployment
  + factor(wording) + factor(chamber) + (1|year),
  family = binomial, data = data) 

model3a_edu <- glmer(eduissue ~ factor(age5)*year20 + n_article
  + factor(gender) + factor(edu3) + factor(citysize) 
  + factor(issue-eduissue+welfare)
  + LDPvote + unemployment
  #+ factor(wording) 
  + factor(chamber) + (1|year),
  family = binomial, data = subset(data, data$year>1972))
#eduissue seems to be not asked in 1972

model4_welfare <- glmer(welfare ~ factor(age5)*n_article + year20
  + factor(gender) + factor(edu3) + factor(citysize) 
  + factor(issue)+ LDPvote + unemployment
  + factor(wording) + factor(chamber) + (1|year),
  family = binomial, data = data) 

#for Table 2
result <- cbind(
  summary(model1_welfare)$coefficients[c(2:5, rep(NA, 7)), 1:2],
  summary(model2_welfare)$coefficients[c(2:7, 18, NA, NA, NA, NA), 1:2],
  summary(model3_welfare)$coefficients[c(2:7, 18, 23:26), 1:2],
  summary(model3a_edu)$coefficients[c(2:7, 18, 21:24), 1:2],
  summary(model4_welfare)$coefficients[c(2:5, 7, 6, 18, 23:26), 1:2]
  )

colnames(result) <- c("nocov_est", "se",
  "w.cov_est", "se",
  "intyear_est", "se",
  "intyear_estEDU", "se",
  "intart_est", "se"
)
rownames(result) <- c("20s", "30s", "50s", "60s",
  "20year", "#article", "unemp%",
  "int20s", "int30s", "int50s", "int60s")
result

#for Figure 2 with model 3
options(digits=2)
result_welfare <- summary(
  model3_welfare)$coefficients[c(2:5, 23:26), 1:2]
  
coef<- c(-0.8, 
         result_welfare[8,1],
         result_welfare[7,1],
         0,
         result_welfare[6,1],
         result_welfare[5,1],
         -0.8,
         result_welfare[4,1],
         result_welfare[3,1],
         0,
         result_welfare[2,1],
         result_welfare[1,1],
         -0.8
)

sd<- c(
  0, result_welfare[8,2],
         result_welfare[7,2],
       0,
         result_welfare[6,2],
         result_welfare[5,2],
       0,
         result_welfare[4,2],
         result_welfare[3,2],
       0,
         result_welfare[2,2],
         result_welfare[1,2],
       0
)

treatment <- c("", "over 60 �~ 20 years", 
               "50s �~ 20 years", 
               "40s �~ 1993: reference",
               "30s �~ 20 years", 
               "20s �~ 20 years", "",
               "over 60", "50s", "40s: reference", "30s", "20s", "")

#install.packages("arm")
library(arm)

par(family="serif")
par(ps = 12)

jpeg("fig2_politicsofaging.jpg", 
  height=960, width=960, res=144)
coefplot(coef, sd, varnames=treatment, 
  xlim=c(-0.6, 0.8), cex.pts=1.6, cex.var=1, main = "", 
  var.las=1, mar=c(1,8,3,3))
dev.off() 


##for Table A3 in the online appendix
options(digits=3)
result02 <- cbind(
  summary(model1_welfare)$coefficients[c(2:5, rep(NA, 7), 1,rep(NA, 14)) , 
    c(1:2, 4)],
  summary(model2_welfare)$coefficients[c(2:7, 18, rep(NA, 4), 1, 8:17, 19:22), 
    c(1:2, 4)],
  summary(model3_welfare)$coefficients[c(2:7, 18, 23:26, 1, 8:17, 19:22),
    c(1:2, 4)],
  summary(model3a_edu)$coefficients[c(2:7, 18, 21:24, 1, 8:17, NA, NA, 19:20),
    c(1:2, 4)],
  summary(model4_welfare)$coefficients[c(2:5, 7:6, 18, 23:26, 1, 8:17, 19:22),
    c(1:2, 4)]
)

rownames(result02) <- c("Age:20s",
  "Age:30s",
  "Age:50s",
  "Age:60s",
  "20yrs:1993 as reference",
  "Media Coverage: Number of Articles",
  "Econ Indicator:Unemployment rate",
  "Age:20s*20yrs",
  "Age:30s*20yrs",
  "Age:50s*20yrs",
  "Age:60s*20yrs",
  "intercept",
  "Gender:Female",
  "Education:High school",
  "Education:Some college or more",
  "Citysize:city w. more than 100k",
  "Citysize:city w. more than 100k",
  "Citysize:Town or village",
  "Main Issues: chose one",
  "Main Issues: chose two",
  "Main Issues: chose all three",
  "Ideology:LDP voting",
  "Wording:Welfare&Elderly Nursing Care",
  "Wording:Medical&Elderly Nursing Care",
  "Election:Upper House",
  "Election:Simultaneious"                       
)

colnames(result02) <- c(
  "nocov_est", "se", "p",
  "w.cov_est", "se", "p",
  "intyear_est", "se", "p",
  "intyear_estEDU", "se", "p",
  "intart_est", "se", "p"
)
result02

#for Table A4 in the online appendix
model5_main <- glmer(issue + welfare ~ factor(age5) + year20 
  + factor(gender) + factor(edu3) + factor(citysize) 
  + factor(wording) + factor(chamber) + (1|year),
  family = poisson, data = subset(data, data$year>1972)
  )

